#!/usr/bin/perl
use strict;
# use warnings;
use Cwd;
use lib getcwd."/bin/bioperl-1.5.2_102";

use Bio::TreeIO;

#####################################################

# my ($USAGE) =  "\n\n\t****** USAGE of preprocess.pl PROGRAM ******\n\n\n\n\tUSAGE: perl $0 <panam.ini file> \n\n\n\n";
# die "$USAGE" if ((scalar(@ARGV))< 1);
# my $option_file = $ARGV[0];
my $option_file = "panam.ini";
chomp($option_file);
my $taxo_group = $ARGV[0];
chomp($taxo_group);

die "\n\n\t=> Cannot find configuration file: $option_file.\n\n" unless (-e $option_file);
die "\n\n\t=> Configuration file, $option_file, appears to be empty!\n\n" if (-z $option_file);
open(OFILE, "<$option_file") || die "Cannot open input file : $option_file\n";

my $usearch;
my $NGS_id_Results;
my $panam_output = "panam_output";
# my $parsing;
my $query_seq; 
my $user_file;
my $preprocess_output = "preprocess_output";
my $seq_F;
my $seq_R;
my $primF;
my $primR;
my $trim =0;
my $path = "bin/uclust3.0.617";
my $us_version = "3.0.617";
####
my $lib = "lib";
my $fast = "$lib/FastTree";
my $align = "$lib/hmmalign";
my $Reference = "Reference" ;
my $build = "$lib/hmmbuild";
####
my $dom;
my %dom;
my $aff_norm;


while (<OFILE>) {
	chomp;
	my $line = $_;

# discard blank line
	if ($line =~ m/^\s*$/g) {
		next;
	}

# discard comment line
	elsif($line =~ m/^\s*#/g) {
		next;
	}
	
# get input options line
        else {
		
		if($line=~m/454_RUN_IDENTIFIER/gi) { 
			$line=~s/454_RUN_IDENTIFIER//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$NGS_id_Results = $line; 
			}
	 		die "\n\tOutput directory for NGS analyses is not defined. Check panam.ini.\n\n" unless ( $NGS_id_Results ne "");
			unless (-d $NGS_id_Results) { mkdir $NGS_id_Results || die "\n\tCould not create Directory $NGS_id_Results $!\n"; }
		}

# 		if (-d $NGS_id_Results."/".$panam_output) { `rm -r $NGS_id_Results/$panam_output` } 
		if (!(-d $NGS_id_Results."/".$panam_output)) { mkdir $NGS_id_Results."/".$panam_output || die "\n\tCould not create Directory $NGS_id_Results/$panam_output $!\n"; }


#		if($line=~m/QUERY_SEQUENCES/gi) {
#			$line=~s/QUERY_SEQUENCES//gi;
#			$line=~s/\t//gi;
#			$line=~s/^\s+//gi;
#			$line=~s/\s+$//gi; 
#			if($line) {
#				$query_seq = $line; 
#			}
#			chomp ($query_seq);
#			# See if the query sequences file exist and contains something
#			die "\n\n\tCannot find query sequence file: $query_seq. Check panam.ini.\n\n" unless (-e $query_seq);
#			die "\n\n\tQuery sequence file, $query_seq, appears to be empty!\n\n" if (-z $query_seq);
#			
#			my $ligne1= `sed -n '1p' $query_seq`;
#			my $ligne2= `sed -n '2p' $query_seq`;	
#			if (($ligne1 !~ /^>/) or ($ligne2 !~ /^[atcguATCGU]/)) {
#				die "\n\n\tYour file $query_seq does not seem to contain sequences in fasta format.\n\n";
#			}		
#		}

		if($line=~m/USER_FILE/gi) {
			$line=~s/USER_FILE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$user_file = $line; 
			}
			die "\n\n\t User file value missed. Check panam.ini.\n\n" unless (($user_file ne ""));
			die "\n\n\t Invalid user file value!\n\n" unless (($user_file eq "own") or ($user_file eq "pooledSamples") or ($user_file eq "eachSample"));	
		}

#		if($line=~m/REFERENCE_BASE/gi) {
#			$line=~s/REFERENCE_BASE//gi;
#			$line=~s/\t//gi;
#			$line=~s/^\s+//gi;
#			$line=~s/\s+$//gi;
#			if($line) {
#				$Reference = $line; 
#			}
#		}

		if($line=~m/FORWARD_PRIMER_NAME/gi) {
			$line=~s/FORWARD_PRIMER_NAME//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$primF = $line; 
			}
		}
	
		if($line=~m/REVERSE_PRIMER_NAME/gi) {
			$line=~s/REVERSE_PRIMER_NAME//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$primR = $line; 
			}
		}

		if($line=~m/FORWARD_PRIMER_SEQUENCE/gi) {
			$line=~s/FORWARD_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seq_F = $line; 
			}
		}
	
		if($line=~m/REVERSE_PRIMER_SEQUENCE/gi) {
			$line=~s/REVERSE_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seq_R = $line;
			}
		}

		############### normalisation
		if($line=~m/NBR_SEQ_NORM/gi) {
			$line=~s/NBR_SEQ_NORM//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$aff_norm = $line; #print "aff --- $aff_norm\n";
			} 
		}
		
		#####################################

		if($line=~m/DOMAIN/gi) {
			$line=~s/DOMAIN//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) { 
				$dom = $line;			
				if (defined $dom ) { #push (@dom, $dom);
					$dom{$dom} = 1
				}
			}
			foreach my $d (keys %dom) {
				die "\n\n\t Domain name is not correct. Check panam.ini.\n\n" unless (($d eq "bacteria") or ($d eq "eukaryota") or ($d eq "archaea"));
			}	
		}
	}
}

die "\n\n\t Invalid number of sequences to pick for the normalization! Check panam.ini.\n\n" if ((defined $aff_norm ) and ($aff_norm !~ /^\d+$/));
if (!(defined $aff_norm)) {$aff_norm = '200' }


print "$Reference\n";

############################################# trimmed profiles path ######################################################
# # # 
if (!(defined $primF) and (defined $primR))  {
	die "\n\n\t Primer forward not defined. Check panam.ini.\n\n";
}

if (defined $primF) {
	die "\n\n\t Primer reverse not defined. Check panam.ini.\n\n" unless (defined $primR);
	die "\n\n\t Primers' sequences missed. Check panam.ini.\n\n" unless ((defined $seq_F) and (defined $seq_R));
	die "\n\n\t Primers' sequences seem to be in the wrong format. Check panam.ini.\n\n" unless (($seq_F =~ /[ATCGUatcgu]/) and ($seq_R =~ /[ATCGUatcgu]/));
	$trim = 1
}

my @tax=();
my $Profiles;

if ($trim == 1) {
	#$Profiles = "Reference/Profiles_".$primF."_".$primR;
	$Profiles = "Reference/Trimmed_Profiles";
	if ((-d $Profiles)) {
		print "Trimmed profiles exist in $Profiles/\n";

	}
}

 
if ($trim == 0) {
	$Profiles= "Reference/Profiles";
}


# 
# # ##############################""
# 

#if (-e $NGS_id_Results."/".$panam_output."/PANAM_Affiliation.txt") {
#		`rm $NGS_id_Results/$panam_output/PANAM_Affiliation.txt`;
#}
#
#`cat $NGS_id_Results/$panam_output/PANAM_Affiliation_*.txt >> $NGS_id_Results/$panam_output/PANAM_Affiliation.txt`;

open (CLADENN, ">".$NGS_id_Results."/".$panam_output."/PANAM_Clades_NN_".$taxo_group.".txt") || die "can not open file";
open (CLADELCA, ">".$NGS_id_Results."/".$panam_output."/PANAM_Clades_LCA_".$taxo_group.".txt") || die "can not open file";

my %sample;
my %taxOTU;
my %leveltax;

########################### clades
if ( -e $NGS_id_Results."/panam_output/Phylogeny/Node_files/file_get_descendent_seq_".$taxo_group){
		open (W, $NGS_id_Results."/panam_output/Phylogeny/Node_files/file_get_descendent_seq_".$taxo_group) ;

		my $i=0;
		my %hh; my %child; my %node_seq; my %boot;
		

################################ /!\ modif : panam_clades par grp taxo 28/3/2012
	if (-e $NGS_id_Results."/".$panam_output."/PANAM_Affiliation_".$taxo_group.".txt") { 
		open (AA, $NGS_id_Results."/".$panam_output."/PANAM_Affiliation_".$taxo_group.".txt");

		my $c; my %hhh; my $bb;  my %bootstrap;
		while (<AA>) {	
			if ($_=~ /^>Read: (.*?)_/) {
				my $aa = $1;
				$sample{$aa}=1;
			}
			if ($_=~ /^>Read: (.*?)\t(.*?)\n/) {
				$c = $1; 
				$bootstrap{$c} = $2;
			}
			else {
				chomp($_);
				$hhh{$c}.= $_;
			}
		}
		close AA; 

		my %clade;my %nbr; my %tax; my %boot;

		foreach my $f (keys %hhh) { 
			if (($hhh{$f} !~ /Metazoa/) and ($hhh{$f} !~ /Other/)) {
				$bb = $bootstrap{$f};
				my @tt = split (/Lowest node:/, $hhh{$f}) ;
				chomp ($tt[0]); chomp ($tt[1]);
	
				################"" %taxOTU pour user_file = pooled et %leveltax pour user_file = eachSample
	
				my @uu = split (/\t/, $tt[1]);
				my $node = $uu[0]; 

				$boot{'node'}{$node} = $bb;
				$clade{'node'}{$node}.=$f.", ";
				$nbr{'node'}{$node}++;
				$tax{'node'}{$node} =  $uu[1];
	
				if ((defined $uu[1])  and ($uu[1] ne "")) {
					$taxOTU{$f}{'LCA'} = $uu[1]; 
					chomp($uu[1]); 

					$uu[1]=~s/\t//gi;
					$uu[1]=~s/^\s+//gi;
					$uu[1]=~s/\s+$//gi;

					my @aqs = split (/;/, $uu[1]);
					my $aqs = $#aqs; 

					if ($aqs == 0) { }
		
					else {
						$leveltax{'LCA'}{$aqs[1]}{$aqs[2]}++
					}
				}
	
		
				my @u = split (/\t/, $tt[0]);
				if ((defined $u[1])  and ($u[1] ne "")) {
					$taxOTU{$f}{'NN'} = $u[1];
					chomp($u[1]); 
			
					my @bqs = split (/;/, $u[1]); 
					my $bqs = $#bqs;

					if ($bqs == 0) { }
					else {
						$leveltax{'NN'}{$bqs[1]}{$bqs[2]}++
					}
				}

###########################################"

				my @t = split (/Nearest neighbor:/, $tt[0]);
		
				my @proc = split (/\t/, $t[1]); 
				my $voisin = $proc[0]; 
	
				$boot{'voisin'}{$voisin} = $bb;
				$tax{'voisin'}{$voisin}=$proc[1]."\n\t\t";
				$clade{'voisin'}{$voisin}.=$f.", ";
				$nbr{'voisin'}{$voisin}++;	
			}
		}

		foreach my $p (keys %{$clade{'voisin'}}) {
			chomp ($tax{'voisin'}{$p}); 
			if ($nbr{'voisin'}{$p} == 1 ) {
				print  CLADENN ">$p\t$tax{'voisin'}{$p}Phylogeny:  ".$taxo_group."_rooted.newick\t".$boot{'voisin'}{$p}."\n\n\t".$nbr{'voisin'}{$p}." neighbor: ".$clade{'voisin'}{$p}."\n\n";					
			}
			else {
				print  CLADENN ">$p\t$tax{'voisin'}{$p}Phylogeny: ".$taxo_group."_rooted.newick\t".$boot{'voisin'}{$p}."\n\n\t".$nbr{'voisin'}{$p}." neighbors: ".$clade{'voisin'}{$p}."\n\n";
			}
		}

		foreach my $p (keys %{$clade{'node'}}) { 
			chomp ($tax{'node'}{$p}); 
			if ($nbr{'node'}{$p} == 1 ) {
				print  CLADELCA ">$p\t$tax{'node'}{$p}\n\t\tPhylogeny:  ".$taxo_group."_rooted.newick\t".$boot{'node'}{$p}."\n\n\t".$nbr{'node'}{$p}." Descendent: ".$clade{'node'}{$p}."\n\n";
			}
			else {
				print  CLADELCA ">$p\t$tax{'node'}{$p}\n\t\tPhylogeny:  ".$taxo_group."_rooted.newick\t".$boot{'node'}{$p}."\n\n\t".$nbr{'node'}{$p}." Descendents : ".$clade{'node'}{$p}."\n\n";
			}
		}
	}
}	
	
#######################################"fin modif : panam_clades par grp taxo 28/3/2012
	
close CLADENN;
close CLADELCA;
